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Abstract. 

Using the non-empirical Variational Induced Breath- 
ing (VIB) model, the thermal properties of periclase 
(MgO) under high pressures and temperatures are in- 
vestigated using molecular dynamics, which includes all 
anharmonic effects. Equations of state for temperatures 
up to 3000K and pressures up to 310 GPa were calcu- 
lated. Bulk modulus, thermal expansivity, Anderson- 
Griineisen parameter, thermal pressure, Griineisen pa- 
rameter and their pressure and temperature dependen- 
cies are studied in order to better understand high pres- 
sure effects on thermal properties. The results agree 
very well with experiments and show that the thermal 
expansivity decreases with pressure up to about 100 
GPa (77=0.73), and is almost pressure and temperature 
independent above this compression. It is also effected 
by anharmonicity at zero pressure and temperatures 
above 2500K. The thermal pressure changes very lit- 
tle with increasing pressures and temperatures, and the 
Griineisen parameter is temperature independent and 
decreases slightly with pressure. 

Introduction 

The density of lower mantle minerals at high pres- 
sures and temperatures is crucial for interpreting seis- 
mological data in terms of mineral constituents. For 
the first time all anharmonic effects up to extreme pres- 
sures and temperatures are being addressed using a non- 
empirical, first principles ionic model. We study MgO 
here since 1) We have an extremely well characterized 
and accurate potential for MgO; 2) it is an end member 
of magnesiowiistite (Mg,Fe)0 thought to be a major 
constituent of the lower mantle; 3) MgO has no phase 
transitions or soft modes in the region of interest and 
therefore should show generic oxide character; and 4) 
MgO is the simplest oxide; other dense minerals whose 
compression is homogeneous should behave similarly. 

Previous studies using non-empirical potentials [Isaak 
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et al. 1990; Wolf and Bukowinski, 1988], empirical po- 
tentials [Agnon and Bukowinski, 1990, Reynard and 
Price, 1990] and thermodynamic models [Anderson et 
al. 1992, 1993] all concluded that the thermal expansiv- 
ity becomes temperature independent at high pressures. 
All studies are consistent with an Anderson-Griiniesen 
parameter 6t = ( oinV ) T ec l ua l to about 5.0 at zero 
pressure which decreases with pressure. How anhar- 
monicity might effect these results and how St behaves 
with pressure and temperature, as well as other thermal 
quantities like Kt, J and Pth is addressed here. 

Method 

We have developed a massively parallel molecular dy- 
namics code for the connection machine CM-5, using 
the VIB model [Wolf and Bukowinski, 1988] and peri- 
odic boundary conditions. We have used supercells of 
64 and 216 atoms, and performed constant E, N and 
V simulations; no scaling was used after equilibration. 
Results of P and T were obtained by averaging over 4- 
10 ps after the system has equilibrated, which took 1-2 
ps. The runs were sufficiently long so that statistical 
errors were smaller than the size of the symbol. The 
Variational Induced Breathing (VIB) potential is based 
on density functional theory which states that all the 
ground state properties are determined by the charge 
density alone, which is modeled here by overlapping 
ionic charge densities. Since oxygen ions are unstable in 
the free state, they are surrounded by a Watson sphere 
of charge +2. The Watson spheres are allowed to re- 
lax spherically at each time step to minimize the total 
energy. In the related PIB model the Watson sphere ra- 
dius is chosen to give the Madelung potential at the ion 
site. The Local Density Approximation (LDA) is used 
to calculate the total energy which is written as the sum 
of the self-energy of each ion, the electrostatic energies, 
the overlap kinetic energy (calculated using the local 
Thomas-Fermi electron gas functional) and exchange 
and correlation energies, which are all functionals of 
the charge density. Further discussion of the method 
can be found in Cohen and Gong [1994], who studied 
the melting of MgO using PIB rather than VIB, and a 
finite cluster rather than periodic boundary conditions. 

Results and Discussion 

We find excellent agreement between our equation 
of state parameters and experimental values (Table 1). 
The temperature variation of the thermal expansivity 
is shown in Fig. la at pressures up to 310 GPa. Since 
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we use classical MD, and quantum thermal and zero 
point corrections are neglected, we expect the results 
to be accurate at temperatures above the Debye tem- 
perature (0£)=945K for MgO). Above the Debye tem- 
perature the thermal expansivity agrees well with ex- 
perimental data. Previous non-empirical studies have 
never included all effects of anharmonicity; here at high 
temperatures and zero pressure we find a nearly linear 
thermal expansivity with temperature compared with 
the diverging form found by Isaak et al. [1990] which 
employed the quasi-harmonic approximation. At higher 
pressures (lOOGPa) and high temperatures we do not 
find significant differences between our results and the 
quasi-harmonic results implying therefore that at high 
pressures anharmonicity does not play a major role. At 
high pressures (Fig. lb) we also find excellent agree- 
ment between the results of this study and experimental 
data. Duffy and Ahrens [1993] measured the thermal 
expansivity at T = 1800K and pressures ranging be- 
tween 160 and 200 GPa to be a constant within exper- 
imental errors. These results agree extremely well with 
our values for the T=1880 K isotherm at these pres- 
sures. We further predict the thermal expansivity to be 
pressure and also temperature independent at higher 
pressures and at higher temperatures as previous non- 
empirical studies have not studied the effects of such 
extreme pressures, up to 310 GPa, on a. The ther- 
mal expansivity decreases with pressure up to about 
100-150GPa (ry=0.73-0.65) at which point it becomes 
pressure and temperature independent. 

Shown in Fig. 2 as a function of pressure is the 
Anderson-Griineisen parameter, St = (dlna / dlnV)T ■ 
For St at ambient conditions, we obtain (5^=4.8, in 
excellent agreement with experiment [Chopelas and 
Boehler, 1989]. With increasing pressure the value of 
the Anderson-Griineisen parameter decreases and be- 
comes temperature independent, much like Agnon and 
Bukowinski's [1990] empirical model. Anderson and 
Isaak [1993] have suggested that St = 6T r) K , with k 
= 1.4 and independent of T and rj. We find that this 
is correct only up to about 77=0.8 (50 GPa), at which 
point this begins to diverge from our results (Fig. 2). 

The quantity <xKt is important in determining other 
thermodynamic properties. The bulk modulus Kt is 
also almost temperature independent at high pressures 
(Fig 3) and is consistent with the above results as 
(§p )t = -j^r ( ) P ■ ^ zero pressure we get very good 
agreement with experiment and at high temperatures 
there is no difference between our results and the quasi- 
harmonic results implying anharmonicity does not effect 
K T . 
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The Gruneisen parameter 7 depends on cxKt via 
7 = , and its temperature and pressure depen- 

dence are shown in Fig 4a and b, respectively. We find 
that the Gruneisen parameter is temperature indepen- 
dent in accord with experimental results and decreases 
slightly as a function of pressure, which is similar to 
other theoretical models [Agnon and Bukowinski, 1990] 
and calculations [Isaak et at, 1990b]. Here too, the ef- 
fects of not including anharmonicity are shown by the 
diverging form of the dashed line corresponding to the 
results of Isaak et al. who used the quasi-harmonic ap- 
proximation. 

The thermal pressure along an isochore is defined as 
P TH = jJ(dP/dT) v dt = Const. + aK T T. We cal- 
culated the pressure dependence of the thermal pres- 
sure at different temperatures above the Debye temper- 
ature. As pressure increases, the thermal pressure in- 
creases, up to about 100 GPa (77=0. 73), at which point 
the thermal pressure becomes pressure independent. In 
order to understand this phenomenon we can look at 
( dl ' a g I p r ^ )t=ol(K t — St). Our results show that while 
K T is slightly smaller than St at zero pressure (K T = 
4.68, (5^=4.8 and K t =St at 20 Gpa), St decreases much 
faster than K T at pressures up to 150 GPa [St decreases 
exponentially and K T decreases linearly). At this point 
St becomes almost pressure independent while K T con- 
tinuously decreases until at 250 GPa K =St ■ 

Conclusions 

These ah initio calculations yield values for the ther- 
mal expansivity which are in very good agreement with 
available experimental data and show that with increas- 
ing pressure the thermal expansivity decreases, up to 
about 120 GPa (ry=0.7) at which pressure it becomes 
pressure and temperature independent. Anharmonicity 
effects are shown to be important only at zero pres- 
sure and temperatures above 2500K. The Anderson- 
Griineisen parameter St follows the same behavior; it 
decreases with increasing pressure and become tempera- 
ture independent at high pressures. The bulk modulus 
is found to be temperature independent at high pres- 
sures and unaffected by anharmonicity. The Gruneisen 
parameter 7 and the thermal pressure are also calcu- 
lated at various pressures and temperatures and are also 
found to be almost temperature and pressure indepen- 
dent, which is consistent with the assumption made by 
Birch [1952] that uKt is independent of density in the 
deep earth. Previous non-empirical studies have never 
found such good agreement with experiment. 
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Figure 1. a). Temperature variation of the thermal ex- 
pansivity at pressures up to 310 GPa. Lines for 0, 100, 
200 and 310 GPa are results from this study, dashed line 
is from Isaak et al. [1990] +'s are data from Touloukian 
et at, and *'s are data from Anderson et al. [1993a] at 
100 GPa. b). Pressure variation of the thermal expan- 
sivity at different temperatures above the Debye tem- 
perature. Lines are results of this study, »'s are data 
from Duffy and Ahrens at 1800K, x's are from Chopelas 
and Boehler [1992] at 2000K. 
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Figure 2. Pressure Variation of St at different tem- 
peratures. The o's are T=940K, +'s are T=1880K and 
*'s are T=2830K. St decreases with pressure and be- 
comes temperature independent at high pressures, x's 
represent Anderson-Isaak proposed fit of St=St r/ 1A . 



Figure 3. The temperature variation of the Bulk mod- 
ulus for different pressures. The +'s are data from An- 
derson et al. [1992], and the solid line is from Isaak et 
al. [1990] using the PIB model and the quasi-harmonic 
approximation. 



Figure 4. a). The temperature dependence of the 
Griineisen parameter at zero pressure. The solid line 
represents results of this study, the dashed line is from 
calculations by Isaak et al. [1990] and the +'s are 
experimental results from Sumino, Isaak and Garvin, 
b). Pressure variation of the Griineisen parameter at 
1000K. Results from this study are shown by o's, +'s 
represent results from calculations by Isaak et al. [1990] 
and *'s are calculations from a model by Anderson et 
al. [1993a]. 



Table 1. Calculated Equation-of-state Parameters for MgO at zero pressure. 





T = 300 K 


This Study 




Exp. 


This Study 


QH PIB a 


QH VIB b 


EM C 


EM d 


940K 


1880K 


2830K 


V [A 3 ] 


18.66 6 


18.70 


18.66 


18.69 


18.73 




19.29 


20.29 


21.38 


K T [GPa] 


160. l f 


153.11 


180.1 


177 


232 


161.3 


132.09 


104.15 


83.93 


r i 
K T 


4.228 


4.68 


4.15 








5.01 


5.44 


5.66 


Krp 




-0.05 


-0.026 








-0.067 


-0.07 


-0.09 


a xlO" 6 [A' -1 ] 


32. h 


38.79 


23.9 




19.4 


31 


44.56 


52.41 


58.75 



a) QH:Quasi-harmonic, Isaak et al. [1990], b) Wolf and Bukowinski [1988], c) EM:Empirical model, Reynard and Price 
[1990], d) Agnon and Bukowinski [1990], e) Mao and Bell [1979], f) Sumino et al. [1983], g) Chang and Barsch [1969], h) 
Touloukian et a/. [1977] 



7 



80 - 



T 60 



> 

§40 

Q. 
X 
<D 

"5 

E 

©20 








^ GPa 



-X- 



i 



.... J&....T.. 



100 GPa. 



200_GPa 



lZ3fo"GPa 



o 



1000 2000 
Temperature [K] 



3000 




J I I I I I I 

100 200 300 

Pressure [GPa] 








100 200 
Pressure [GPa] 



300 



1000 



800 



o 

Q_ 
O 

l -|500 

(0 

"D 
O 

^400 



■3 
m 



200 z 



r- 1 


1 1 

310 


1 

CPn 
\J r u 




.Q. 


_e 




200 


GPa 











100 


GPa 








PIB-n 






^EXP 


GPa J 

i 


1 


1000 


2000 


3000 


Temperature [K] 





i r 



P = GPa - 



+ + + + + + 

jK-^.^ <j,_ »R — 



MD 



-o 



PIB 



j i_ 



o 



1000 2000 
T [K] 



3000 




1 00 200 

Pressure [GPa] 



